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Abstract 

This paper develops methodology for local sensitivity analysis based on directional deriva¬ 
tives associated with spatial processes. Formal gradient analysis for spatial processes was 
elaborated in previous papers, focusing on distribution theory for directional derivatives 
associated with a response variable assumed to follow a Gaussian process model. In the 
current work, these ideas are extended to additionally accommodate a continuous covariate 
whose directional derivatives are also of interest and to relate the behavior of the directional 
derivatives of the response surface to those of the covariate surface. It is of interest to assess 
whether, in some sense, the gradients of the response follow those of the explanatory variable. 
The joint Gaussian structure of all variables, including the directional derivatives, allows for 
explicit distribution theory and, hence, kriging across the spatial region using multivariate 
normal theory. Working within a Bayesian hierarchical modeling framework, posterior sam¬ 
ples enable all gradient analysis to occur post model fitting. As a proof of concept, we show 
how our methodology can be applied to a standard geostatistical modeling setting using a 
simulation example. For a real data illustration, we work with point pattern data, deferring 
our gradient analysis to the intensity surface, adopting a log-Gaussian Gox process model. 
In particular, we relate elevation data to point patterns associated with several tree species 
in Duke Forest. 

Keywords: Gauchy Process, Directional Derivative, Gaussian Process, log-Gaussian Gox 
Process, Matern Gorrelation Function 


1. Introduction 

Increasingly, data is being collected at geo-referenced locations. For a region of interest 
D, the set of conceptual responses {V(s) : s G D} can be viewed as a realization of a 
random surface, observed at a finite set of locations. While covariate information may 
explain a substantial portion of the variation in response, there is often underlying spatial 
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structure that is difficult to measure. Inference on this spatial structure can be made via 
the parameters in a spatial process model. Under these models, prediction at unobserved 
locations, or kriging, is possible, enabling interpolation across the region. 

Spatial regression models commonly assume a linear relationship and make inference 
based on the coefficient assigned to the covariate. This coefficient describes the expected 
change in response given a unit change in covariate, thus providing a global measure for the 
sensitivity of the response to the covariate. However, it is expected that the relationship 
between the variables may vary locally over the study region. Such local, or second-order, 
behavior can be studied through spatial sensitivity or gradient analysis. 

A spatial gradient analysis will enable spatial examination of a response variable’s sen¬ 
sitivity to a covariate across the region of interest. The sensitivity of the response variable 
may vary based on the rate of change for the covariate or due to additional unaccounted 
for factors, resulting in areas where the relation ship appears weaker or stronger. Mode ls 


allow ing for spatially varying coefficients, such as fiFotheringham et al.l. 120021: iGelfand et ah 


2003|), provide some similar inference but assume a more complex model structure. The 


methodology proposed here assumes a standard spatial linear regression model but provides 
a post model htting framework for examining the variation in the response’s sensitivity to 
the covariate. 

In ecology such sensitivities are typically discussed when relating plant characteristics 
to climate. For instance, researchers are increasingly interested in characterizing how abun¬ 
dance and frequency of tree species relate to changes in variables such as temperature and 
precipitati on, in order to l earn about the expected effects of climate change o n range distribu¬ 
tions (e.g. Thniller et ah . 2004; Can ham and Thomasl . 2010 : Thomas, 2010 1. These analyses 
focus on a comparison of sensitivities across species, while little consideration is given to how 
sensitivities may vary spatially within a given species. Our approach is aimed at the latter 
question. 

The interpretation of the coefficient in a spatial regression as a global gradient, dE{Y (s))/ 
dX{s), inspires consideration of local sensitivities through directional derivative processes, or 


spatial gradients. Spatial gradients under Gaussian processes were elaborated bv iBaneriee et ah 


(l2003bl l to address the rate of change of a spatial surface at a given point in a given direction. 
Their paper dehnes directional derivative processes with corresponding distribution theory 
to enable interpolation across a region. The gradient distributions are fully determined by 
the spatial model parameters, allowing all gradient analysis to occur post model htting. Dis¬ 
tributions for derivatives of Gaussian processe s have also been d i scusse d in the context of 


observed derivatives of fun ctions (e.g. OHaganl 19921: Solak et ah . 2003 1. as well as for ran¬ 


dom helds more generally ( Adlerl, 1981). In all of the previous work with Gauss i an pr ocess 


directional derivatives (e.g. Baneriee and Gelfandl. 2003 . 2006 : Maiumdar et ah . 2006), the 


researchers have considered the rates of change of a response surface with the mean surface 
modeled as a linear function of a set of hxed covariates. In contrast, to accommodate the 
desired spatial sensitivity analysis, we will assume a single covariate of interest whose surface 
is spatially smooth such that it too can be treated as a realization of a stochastic process. 
The behavior of the response and covariate processes, as well as their associated derivative 
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processes, are then considered jointly and functions of these derivatives can be explored. 

The contribution of this paper is to extend the existing spatial gradient theory to accom¬ 
modate spatial sensitivity analysis by modeling the response and covariate jointly. Working 
within a hierarchical Bayesian modeling framework, corresponding gradients for the spatial 
surfaces can be sampled simultaneously from the joint predictive distribution post model £t- 
ting. Under a significant regression relationship it is not sensible to investigate the gradient 
behavior of the surfaces marginally. Suitable comparison between the gradient surfaces illus¬ 
trates how sensitive the response surface is to the covariate surface, as well as the strength 
of this relationship. The former is accomplished through comparison between the directions 
of the maximum gradient at a given location; the latter requires consideration of their di¬ 
rectional derivatives relative to one another. In particular, we introduce two new spatial 
processes, a local directional sensitivity process and a spatial angular discrepancy process. 
These inferential tools are devel oped and carried out on simulated data in the c ontext of a 
customary geostatistical model ( Baueriee et ah . 2003al: Cressie and Wikle . 2011 1 as well as 
with an ecological dataset where we connect point patterns of trees with elevation. 

In Section [2] the formal distribution theory for the spatial gradients is extended to the 
multivariate case. Section E] outlines the modeling framework for our examples. This section 
also defines the two processes of interest, namely the local directional sensitivity process 
and the spatial angular discrepancy process. Section 0] provides a simulated example with a 
multivariate Gaussian process setup as a proof of concept. Section [5] provides an analysis of 
point pattern data from Duke Forest, extending the analysis techniques to a non-Gaussian 
response; the intensity of the point pattern, modeled using a log-Gaussian Gox process, 
is explored by employing a spatial gradient chain rule. Finally, Section [6] summarizes the 
contributions of the paper and suggests future work. 


2. Distribution Development 


In t his section we review the dehnitions and distributions presented in iBaneriee et ah 


fl2003bl) and extend these ideas to consider a multivariate Gaussian process. We assume lo¬ 


cations s G 2-dimensional Euclidean space, however extension to a generic d-dimensional 
setting is straightforward. The process is assumed, for convenience, to be (weakly) stationary 
such that the covariance function, U(s')), depends only on the separation vector 

(5 = s — s'. In fact, in our examples we adopt isotropic covariance functions that depend 
only on the length of the separation vector, ||5||. 

Gonsider two surfaces {(y'(s),X(s)) : s G drawn from a joint Gaussian process 
specihed such that W(s) has constant mean, say Oq, and covariance function G{S). Given 
df(s), U(s) has mean /3X(s) and covariance function K{d). Observed at a set of locations 
Y = (Y(si),..., Y(s„)), we write: Y|X ~ iV(/3X, iF(-)) with X ~ Y(q;o, G'(-)). Gonsidered 
jointly, we have: 


Y 

X 


N 


ao/51 

Uol 


(5G{-) 


^G{-) 

G(-) 
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where G{-) and K{-) are matrices of the covariance functions with entry i,j evaluated at 
S = Si — Sj. 

We follow the notation and theory in iBaneriee et all fj2003bl) . Suppose mean square 
differentiable processes f^(s) and W(s). That is, for T(s), at Sq there exists a vector Vy(so) 
such that for any scalar h and any unit vector u, Y (so + hu) = Y (sq) + h.u^Vy(so) +r(so, hu) 
where r(so, hu) —)■ 0 in the L 2 sense as h —)■ 0. Similarly, for W(s). 

In particular, dehne the hnite difference processes at scale h in direction u: 

y(s + hu) — Y (s) 


I'ua(s) = 
X„.h(s) = 


h 

df (s + hu) — X(s) 


where u is a unit vector. Taking the limit as h tends to 0 , IBaneriee et al.l fj2003bl) dehne the 
directional derivative processes in the direction u: 

D„y(s) = hmy,,,(s) = u'Vy(s) 
h—>-0 

BuX(s) = limXu^h(s) = u'Vx(s) 

h—>-0 

where Vx(s) = (Dej^X(s), Dg 2 X(s))' is the vector of directional derivatives in the orthonor¬ 
mal basis directions ei = (1, 0) and 62 = (0,1) for We can study the directional derivative 
processes for any u by working w ith the basis set Vy(s) and Vx(s). 

From Banerjee et ahl f 20031)1 ) . we know that if f^(s) and X(s) are stationary Gaussian 
processes, then the resulting marginal distributions involving the directional derivatives will 
be stationary Gaussian processes as well. Note, isotropy in the Y (s) process does not induce 
isotropy in the B^Y (s) process; only stationarity will be inherited. Similar to the discussion 
in their paper, we know by linearity that (y(s), X(s), yu^/i(s), Xu,h(s))' will be a stationary 
multivariate Gaussian process. And then, by a standard limiting moment generating function 
argument, (y(s), A(s), T)uy(s), Zi)uX(s))' will also be a stationary multivariate Gaussian 
process. 

To explicitly provide the joint distribution, we derive the cross covariance structure by 
examining pair-wise covariances between the response and covariate processes and their 
directional derivatives. For notational convenience, write the marginal covariance function 
of T(s) to be K(-) = K(-) Y /3^G(-). Assume the F'(s) and X(s) processes are mean zero, 
setting tto = 0 , since in practice the gradients are calculated for the mean zero residual 
process. If F^(T(s)) = 0, then B(BuY(s)) = 0, so the joint processes will all be mean zero. 
We calculate the covariances associated with the directional derivatives by taking the limits 
of the covariances corresponding to the analogous hnite di fference process . _ 

The covariances for the response surface are derived in Baneriee et ahl f 2003bl) : 


c'o^;(yu,/^(s),yu,/^(s')) = 


2K(S) - K(S + hu) - K(S - hu) 


Gov{B^Y{s),B^Y{s')) = lim Gon(W,ft(s), Fu,h(s')) = -u'Bj^u 

h^O 
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C'ot;(F(s),F„,;,(s')) 

Cov{Y{s),D^Y{s’)) 


K{6 - hu) - K{S) 
h 

lim Cot;(y(s), Yu,,,(s')) 


DuK{-S) 


and covariances for the covariate surface are analogous: 


CoviX^,,,is),X^4s')) 

Cov{D^X{s),D^X{s')) 


2G{S) - G{S + hu) - G{d - hu) 

J? 

liniCon(Xu,,j(s),Xu,,j(s')) = -u'ficu 


C'on(X(s),X„,;,(s')) 


- hu) - ^(5) 
h 


C'on(X(s),D„X(s')) = C'on(X(s),X„,,(s')) 


DuG{-d) 


To fully describe the joint distribution we derive the covariances between response and 
covariate surfaces similarly: 


C'on(F(s),Xu,;,(s')) 

C'on(F(s),DuX(s')) 


/3G{S - hu) - /3G{d) 
h 

limC'on(y(s),X„,;,(s')) = f3D^G{-S) 


Gov{X{s),Y^,h{s')) 

Gov{X{s),D^Y{s')) 


/3G{d - hu) - /3G{d) 
h 

limCon(X(s),Fu,fe(s')) = l3DuG{-S) 


Gov{XuA^),Y^,,,{s') 

Gov{D^X{s),D^Y{s’)) 


2/3G{S) - /3G{d + hu) - /3G{d - hu) 


h2 

limCon(Xu,ft(s),yu,,j(s') = -Pu'Qg{S)u 


where = d‘^G{6)/d6id6j and DuG{d) = limft^o(G(^ — hu) — G{6))/h. 

Relationships between the response surface, the covariate surface and their correspond¬ 
ing directional derivative surfaces can be described through the 6-dimensional multivariate 
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stationary Gaussian process Z(s) = (y(s),X(s) , Vy(s), Vx( s))'. Using the covariances cal¬ 
culated above, the associated cross-covariance matrix for Z will be: 


Vz(S) = 


/ 

k{d) /3G{d) 

-vk{dy 

-/3VG((5)' \ 




/5G((5) G{6) 

-(3VG{6y 





Vk{d) I3VG{6) 

-Hk{S) 

-PHg{5) 



V /3VG((5) VG(5) 

-/3Hg{6) 

-hg{s) y 



/ 

K{d)+kG{6) 

(5G{5) - 

-ViF((5)'-/52VG((5)' 

-/3VG{Sy \ 


^G{5) 

G{6) 

-/5VG(5)' 

-VG{6y 


l> 

+ 

l> 

> 

-Hk{5)-(3^Hg{S) 

-I3Hg{5) 

V 

/5VG(5) 

VG(5) 

-/3HG{d) 

-hg{s) y 


where VK(d) is a 2 x 1 gradient vector associated with K{d), and Hk{S) is the 2x2 Hessian 
matrix associated with K{6). 

For (5 = 0, we have a block diagonal local covariance matrix: 


Uz(0) = 


f K{0)+/3^G{0) /3G(0) 
/3G(0) G(0) 

0 0 
\ 0 0 


0 ' 0 ' \ 

0 ' 0 ' 

-Hk{Q)-I3‘^Hg{Q) -/5ifG(0) 


Thus, at a location s, the directional derivative surfaces will be correlated with one another, 
but neither will be correlated with either of the data surfaces. Intuitively, this makes sense 
since we would not expect the level of the surface at a location to be correlated with the rate 
of change at that location. Of course, since (X(s),y(s))' is a bivariate Gaussian process, 
the correlation between the rates of changes is not surprising. 

The Matern covariance is adopted below. It depends on a smoothness pararneter v which 
directly controls the mean square differentiability of process realizations fjSteinl . Il999[) . This 
is convenient since, again, the U(s) and X(s) processes must be mean square differentiable 
for their associated directional derivative processes to be well dehned. If we let K{-) and 
G{-) be Matern with u > 1 then they are once (but not twice) mean square differentiable, 
and, if 1 / = 3/2, the covariance functions are of the closed form + (/)| |(5| |) exp(—(/)| |<5| |). 

We denote the parameters of K{-) as cr^ and and the parameters of G(-) as and <px- 


y 'f'y 

Under the Matern covariance the components of the cross-covariance matrix will be 

yexp{-(py\\6\\){l-(j)y5f/\\6\\), (Ffx(<5))p 
Then, we have for (5 = 0: 

\ 


ViF((5) = -a2(/)2exp(-(/)j^||(5||)(5, {HK{S))ii = 

= (^y<Plexp{-(py\\S\\)5i6j/\\6 

/ 


and similar for G{- 


Vz{0) = 




(3a\ 


0 

0 


0 

0 


{<ryy + 


( 1 ) 


0 ' 0 ' 

0 ' 0 ' 

Vo 0 {I 3 al 4 >l)l 2 J 

where G is the 2x2 identity matrix. As above, De^Y{s) and D^^Xls) will be correlated 
with one another, and similarly De 2 Y{s) and De^Xls) will be correlated with one another, 
but all other pairings of the directional derivatives will be uncorrelated. 
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3. Model Fitting and Inference 

3.1. Sampling Method 

Following the modeling of the previous section, let K{5) = ayPy{6) and G{d) = a^p^iS), 
where the px and Py are valid two-dimensional correlation functions. We work with the 
Matern class of covariance functions parameterized by (j) and u with u > 1. 

Let 6 = {ao, Po, l3i,al,ay,(j)x,(j)y,i'x,i^y)- For locations Si,...,Sn, the overall likelihood 
can be written in terms of the conditional likelihoods 


L{0;Y,X) oc L{Y\e,X)L{X\0) 


X exp 

X exp 


2<7j 

1 


(X — aol)'Rx^{4>x, i^x)(X — aol) 

(Y - (/3ol + fliX)yR-\(j)y, Uy)(Y - (/3ol + /5iX)) 


where Y , Y(sn)) , (^Rxl^cj^x^ ^x^ij Px{,^i Sjj 0 X 5 and (^Ry(^(l)y, i^y^ij 

Pj^(sj — Sj] (j)y,Uy). The likelihood could be equivalently written in its joint form, but the 
conditional form is more conducive to interpreting and implementing the gradient analysis. 

We see that we have a low dimensional parametric model, with 0 only 9 dimensional. We 
utilize fairly non-informative priors for its components. For example, vague normal priors on 
(tto,005 01)5 vague inverse Gamma priors on (cr^ 5 cr ^)5 vague Ga mma priors on {(f)x,4>y), and 
U ( 1 , 2 ) priors on (z/a,, Uy). The prior on u follows the suggestion of ISteinI (119991 ) and others who 
observe that distinguishing v = 2 from v > 2 would be very difficult in practice. This model 
is straight forward to ht in its conditional form, for example using the ‘ spBayes ’ package in 
R ( Finley et al.l . 20071) . Thus, assume we now have posterior samples 00 Z = 1 ,..., L, from 


/(^|Y,X). 

Once we have posterior samples of the parameters, we draw samples of the gradient 
vectors using composition since the posterior predictive distribution /(Vy,Vx|Y,X) = 
//(Vy,Vx|Y, X,6)f{6\Y,X)d6. The cross-covariance matrix derived earlier allows us to 
immediately write the joint multivariate normal distribution given 0, which can be eval¬ 
uated at each sample 00 Based on this joint distribution, standard multivariate normal 
theory allows us to write down the desired conditional distributions needed to draw from 
the predictive distribution. 

For an unobserved location Sq, obtaining draws of the gradient vectors is again done via 
the predictive distribution. The cross covariance matrix derived enables us to write the joint 
distribution, from which we can derive the conditional distribution /(Vy (sq), Vx(so) | Y, X, 0) 
If interest is also in the values of the Y (s) and X (s) surfaces at the new location, we would 
derive the conditional distribution /(Y(so), X(so), Vy(so), Vx(so)|Y, X, 0), which is again 
straight forward given the cross covariance matrix and allows joint prediction of the surfaces 
and their gradients at the new location. 

If we want Y (s) and X(s) to be adjusted based on some hxed covariates, then we simply 
introduce such covariates into the mean functions of the model. We create a spatial random 
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effects model: 


^ — /^O + + '^y{^yiy + + e(s) 

X(s) = 00 + T^(s)' 7^ + w^(s) 

where Ty(s) and Ta;(s) are vectors of covariates used to explain the Y (s) and X(s) surfaces 
respectively, with coefficients jy and 'y^; Wy{s) and Wx{s) are independent mean-zero station¬ 
ary Gaussian processes with parameters al,ay,(j)x,(t>y,i^x, as before; and e(s) is a Gaussian 
white-noise process with variance intended to capture measurement error or microscale 
variability in the response. The X(s) process is assumed to be a fully spatial model (no 
nugget effect), such as might be used for elevation, temperature, or pollutant level. 

With 6 now extended to 0 = (ao,/do/^i) Ti) WxWy) the likelihood for 

locations Si,..., Sn above can be trivially revised. 

Prior selection for the parameters will be similar to the previous example. Again, this 
model can be implemented using ' spBayes', and draws of the gradients will rely on the 
posterior predictive distribution, which can be calculated as before using the derived cross¬ 
covariance matrix. 


3.2. Local Directional Sensitivity Process 

At a given location s there may additionally be interest in the ratio of directional deriva¬ 
tives, 

DuY{s)/DuX{s), corresponding to the relative rates of change in the two surfaces in di¬ 
rection u. This quantity is analogous to dy/dx in more standard calculus applicat i ons a s 
well as sensitivity functions studied in sensitivity analysis (jTomovic and Vukobrato^ . 1972 ). 
With this in mind, we refer to the resulting spatial process as the local directional sensitivity 
process. The choice of u will depend on the application being considered. For example, this 
direction may correspond to latitudinal direction, an elevation direction, or to an environ¬ 
mental feature expected to impact the response. Large values of this process would suggest 
that the change in covariate surface has a high impact on the the response surface. 

The joint process defined in Section [2] can be equivalently represented as a spatial random 
effects model (excluding nugget effects): 


F(s)|X(s) = l3o + f3iX{s) +Wy{s) 
X(s) = ao + Wx(s) 

Wy(s) ~ gp(o,a:(5)) 

n;,(s) ~GP(0,G((5)) 


where Wy(s) and Wx(s) are independent processes. With this notation we can write the 
unconditional response surface and corresponding directional derivative process as follows: 


V (s) — /3o + /5l«0 + A'^x(s) + UJy(s) 

DuV(s) = /3iDuWx(s) + DuWy(s) 
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The local directional sensitivity process can then be written 


DuY{s) 

DuX{s) 


- = Pi + 


DuW^ 


DuWc, 


(2) 


We see that the multiplicative parameter, /3i, dehning the overall relationship between the 
X(s) and T'(s) processes serves to center the local directional sensitivity process. 

As mentioned in Majumdar et ah 2006, if we consider Y (s) = /5o ++e(s), then one 
could write (3i = dE(Y{s))/dX{s); i.e., /3i is describing the rate of change in E{Y{s)) relative 
to changes in df (s). Again, at any location the directional derivative ratio will be centered at 
the global (non-directional) derivative ratio dE(Y{s))/dX{s) plus some (directional) spatial 
noise. In this way, the directional derivative ratio process is describing the spatial variation 
in the relative rates of change between X( s) and F (s). This is analogous to modeling 
adopting spatially varying coefficients, /9(s), f Gelfand et ah . 2003|l but is arguably a simpler 
context since the derivatives require no additional model fitting. In addition, consideration 
of directional perspectives using this gradient approach allows for inference distinct from 
what one can learn from a non-directional /3{s) parameter. 

By noting that D^iWy{s)/D^j^Wx{s) is a ratio of independent mean zero normal random 
variables, at each location the directional spatial noise is a Cauchy random variable with 
scale equal to the ratio of the respective standard deviations, SD{DuWy{s))/SD{DuWx{s)). 
If u = (mi, M 2), then SD{D-nWy{s)) = ifj^(0)i 1) -|- { 0 ) 2 , 2 )- If Y (s) is isotropic, 

then Hk{ 0) = Coh ( Baneriee et ah . 2003bh and SD(DxiWy{s)) = \/u'l + U 2 \/—Co = y'—C q. 
If K{-) and G{-) are Matern with u = 3/2, then SD{D^Wi{s)) = atcpt, and the scale for 
the Cauchy distribution will be (Jy(j)y/o'x((>x- When the respective standard deviations are 
equal, the scale will be 1 and the directional derivative ratio will have a standard Cauchy 
distribution. 

In fact, the collection of directional derivative ratios form a well defined spatial stochastic 
process which would naturally be called a spatial Cauchy process (see Appendix A for de¬ 
tails). For any set of locations Si,..., the joint distribution is well defined. For example, 
consider two locations s and s', simplifying the notation for clarity: 




D^jWy (s) 

Eu^xi^ J 


DuWy{s') ^ ^ ^ 

< c, ^ < r2) = p {— < ri, — < ra) 


DuWx[S') ■ mi 1712 

= P{ni < rimi, n 2 < 721712, mi > 0 , m2 > 0 ) 

-|- P{ni < TilTli, 172 > 721712, Tfli >0, m 2 < 0) 

-I- P{ni > rimi, M2 < 72m2, mi <0, m2 > 0) 

-I- P{ni > 7imi, 172 > 721712, mi < 0, 1772 < 0) 

In turn, each of these terms can be computed as an integral involving normal densities. For 
example, the first term can be written as follows: 


F(ni < 7imi,172 < 721712,7171 > 0,1712 > 0) = 
roo roo pr\m\ pr2fn2 

/ / / / fK{ni, n2)fG{mi, m2)dn2dnidm2dmi 

' 0 Jo J —oo J —oo 


(3) 
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where /ii:(ni,n 2 ) = fK{DuWy{s),DuWy{s')) is the bivariate normal density for the DuWy{s) 
Gaussian process with parent covariance function K{-) evaluated at s and s', and similar 
for Re call that a univariate Cauch y distribution can be dehned as a scale 

mixture of normals ( Andrews and Mallowsl . 1974I) . If we write the corresponding CDF as 
nr) = ri:: ^(j){n)(j){m)dndm, then the form in Equation [3] is evidently similar and can 
be regarded as a bivariate analogue. In this way, the spatial Cauchy process is dehned such 
that at any location the distribution is a univariate Cauchy distribution and for any set of 
locations the distribution is a sum of integrals of a similar form. 

In the data analysis examples we consider a hxed direction u and draw samples of the 
directional derivative ratio at each location across a region. Although we cannot show mean 
square continuity for this surface, by imposing additional smoothness conditions on the co- 
variance functions for Wy{s) and w-ris), w e can argue that this surface will be almost surely 


U x\^l>i VV’ 

contin uous using results from lKentl fllQSOl ) . following the development in iBanerjee and Celfand 
( 2003 1. We omit the details. 


3.3. Spatial Angular Discrepancy Process 

At any given location there may be interest not only in the magnitude of the gradients in 
various directions, but also in the direction at which the maximum gradients are achieved. 
Stronger alignment between the directions of maximum gradients would suggest a stronger 
relationship between the response and covariate surfaces. 

Consider the covariate process W(s). At location s the maximum gradient will be 
achieved in the direction described by the unit vector = Vx(s)/||Vx(s)||, and the di- 
rectional deriva t ive in the direction of maximum gradient will be X{s) = ||Vx(s)|| 


(iBaneriee et ahl. l2003b|l . We can similarly consider the direction of maximum gradient for 
the response surface which will occur for Uy = Vy(s)/||Vy(s)||. In applications, researchers 


may be interested in the behavior of the response surface, or its relative rate of change, in 
the direction u^. That is, there may be interest in D^* Y(s) = Vjv:(s)'Vy(s)/||Vx(s)|| as 
well as D^*^Y{s)/D^*^X{s) = (Vx(s)'Vy(s))/(Vx(s)'vJ(s)). 

At a location s the magnitude of the maximum gradient, ||Vx(s)||, will be the square 
root of a sum of squared independent normal random variables; as such, this quantity will 
have a Chi distribution with d = 2 degrees of freedom, possibly scaled by some factor. If the 
process is isotropic, then ||Vx(s)|| will have a Chi distribution with d = 2 degrees of freedom 
scaled by SD{Dej^X{s)) = SD{De 2 X{s)). For the Matern covariance {u = 3/2) structure 
this scaling factor will be equal to ax(l)x- 

The unit vector describing the direction of max gradient for the covariate surface can 
equivalently be described by an angle 6x{s) such that tan(6'x(s)) = D(oq)X(s)/D(i o)X(s), 
and similarly for the response surface. The angle 9 can take values from —tt to vr, so inversion 
of tan must be done with care. Since tan has a period of only tt, the inverse function is 
typically taken to be arctan* 


^arctan*(S'/C) is arctan(5'/C') if C > 0, S' > 0; 7r/2 if C = 0, S > 0; arctanlS/C) -|- tt if C < 0; 
arctan(S/C') -I- 27r if C > 0, S < 0; and undefined if C = 0, S = 0 ( .Tammalamadaka and Senguntal . 2001 ). 
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As smooth functions of a well defined spatial process, the directions of maximum gradient 
(6*x(s), 6*y(s)) dehne a bivariate projected Gauss i an pr ocess, analogous to the univariate 
projected Gaussian processes described in Wand (2013 ). Marginally dx(^) and 6'y(s) will 
each be a projected Gaussian process f Wand. 20131) . Assuming the Matern (z/ = 3/2) 
covariance structure we can derive the joint distribution of the two angles at a given location 


s: 


C 


A^(s)0(O) | 


f{9xis),9y{s)) = 


where C = 


ac{ac—A^{s)) {ac—A^{s))^/^ 
A^(s)</)(0) I y/^A{s) 


L(0,0, 


A^is) 


) + « . 


i(s) > 0 


C'(^1A£2|^ + _v|L^( 0.5-L(0,0,G^))+^), A(s)<0 

\ ac(ac—A^{s)) {ac—A^(s))^/^ \ ^ ^ ^ \ ac j ac j ^ ^ ' 


m 

ac 


1^1 = Al&Ayyf. a = i/ialAy), c = 

A(s) = cos(6*x(s) — 9y(s)), and L(0,0,p) is the zero mean bivariate normal cdf with 
correlation p and standard deviations equal to 1 evaluated at (0,0). A brief derivation is 
provided in Appendix B. In addition, it is straightforward to show that at each location s 
the angles, 9x{s) and 0y(s), will marginally be uniform on (—7r,7r). 

The above bivariate density is plotted in Figure [T] for /3 = (±0.05, ±0.5, ±1)', holding 
all other parameters constant: Cy = 1 = (t^ and (j)y = 1.05 = (p^- When /3 > 0 the mass 
is concentrated around {9x,9y) pairs that are equal; when (3 < 0 the mass is concentrated 
around pairs where 9x = 9 y — tt. When (3 = ±0.05 the relationship between X(s) and ±(s) 
is weak, and the density is roughly uniform over all angle pairs. As the magnitude of /3 
increases, the mass becomes increasingly concentrated around these respective values. 

To compare the directions of maximum gradient calculated from the posterior distribu¬ 
tion, we compute a “discrepancy” for 6*x(s) and 6*y(s). Dehne: disc{s) = 1 — cos(6'x(s) — 
9y{s)). As a smooth function of a spatial process, this discrepancy is a well dehned spatial 
process which we refer to as the spatial angular discrepancy process. When the maximum 
gradients occur in identical directions this process will have a value of zero, when they occur 
in opposite directions the process will have a value of two. In our analyses we consider the 
process across a region, plotting the posterior median surface. 

In some areas the X (s) or Y (s) surface may be quite hat and there will be no direction 
with a gradient magnitude substantially larger than in the other directions. In these areas 
a large angular discrepancy may not be as meaningful as it would be in areas with larger 
gradient magnitudes. For this reason it is useful to examine these plots in tandem with plots 
of the local directional sensitivity process described in Section 13.21 

3 . 4 . Extensions to Non-Gaussian Data Models Using the Chain Rule 

One can imagine modeling scenarios where the response will be non-Gaussian, for in¬ 
stance, binary or zero-inhated data. In these cases, a latent Gaussian surface may be uti¬ 
lized and spatial gradients can still be considered. The resulting gradient surfaces will 
be for the latent Gaussian process, not necessarily the response of interest; however, dif¬ 
ferentiable transformations will allow for i nference on the non-Ga ussian response surface. 


Recall the chain rule result presented in iMaiumdar et al.l (120061) : for g{-) diherentiable 
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| 3 = 0.05 


| 3 = 0.5 



Figure 1: Bivariate density for (0x(s), 0v(s)) at a fixed s for varying values of /3. All other parameters are 
set to the values used for simulation in Section H) 


on and iy(s) = 5 f(l^(s)), the directional derivative Ziluhh(s) exists and is given by 
-Duhh(s) = Dug{V{s)) = g'{V{s))D^iV{s). Here H(s) is the latent Ganssian snrface and 
the fnnction g{-) describes its relationship to the non-Ganssian response of interest. 

When analyzing the Dnke Forest data in Section O we work with a log-Gaussian Gox 
process model for the intensity snrface. To stndy gradients associated with the intensity 
snrface itself, we can nse the fact that, if H(s) is onr response and Z(s) = logy(s), then 
D,y(s) =exp(Z(s))Du^(s). 

In the case of binary data we can write a model with response H(s) = 1 if Z(s) > 0 and 
H(s) = 0 if Z{s) < 0, where Z{s) is a latent Ganssian process centered at a linear fnnction of 
a spatial covariate Wfs). We have a spatial probit regression since P(Y{s) = 1) = <h(Z(s)) 
( Heagertv and Lele . 19981) . There may be interest in identifying areas where the transition 
between regions of high absence probabilities to regions of low presence probabilities is rapid 
(or vice versa). Gonsideration of the local directional sensitivity process associated with the 
probability snrface P(s) = P{Y{s) = 1) and the covariate surface W(s) would provide an 
avenue to answer these kinds of questions. Now, with the function g{-) = $(•), we obtain 
D^P{s) = <P{Z{s))D^Z{s). 


4 . Simulation Example 

We consider a simulation example to explore the behavior of gradient quantities in a 
controlled setting. From the foregoing development, in the context of spatial gradient analysis 
the quantities of interest will be directional derivatives. These derivatives are unobservable 
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even in a simulation study with known parameters. Thus, assessing inference performance 
with regard to these quantities requires some novelty. 

Recall that the gradient processes describe the shape and behavior of spatial surfaces. 
With a simulated dataset, we are able to draw a realization of the Gaussian process over a 
larger number of locations, allowing for fairly detailed understanding of the spatial surfaces. 
To avoid any unfair advantage that might come from the increased sample size, we use only 
a subset of the locations to £t the model and reserve the full set of locations for assessing the 
quality of our inference. Contour lines highlighting the shape of the surface are interpolated 
using the full set of locations. Conclusions made using gradients are then compared to those 
suggested by the contour lines to examine performance. 

4 . 1 . Data and Model 

We simulate ^(s) a realization from a mean zero Gaussian process on [0,10] x [0,10], 
and R(s)|X(s) a realization from a Gaussian process with mean /5X(s). We use Matern 
covariance functions setting u = 3/2 in order to capitalize on the resultant closed form. 
We simulate the Gaussian processes assuming Matern covariance functions with parameters 
(f)x = 1.05 = CT^ = 1 = (j^, /So = 0 = ao, and (di = 0.5. We draw a larger realization at 
2000 locations, to allow for finer knowledge of the underlying surface, from which we consider 
a subset of 200 locations to be our “observations”. 

Treating {Y (s), X(s))' as a multivariate Gaussian process, we £t a coregionalization model 
using the conditional parameterization, as described in Section [2l The model is htted using 
the 'spbayes’ package in R by hrst htting parameters for W(s), then htting parameters 
for y(s)|W(s). We obtain 2000 samples after a burn-in of 500 iterations and a thinning of 
every fifth iterate. We assume v = 3/2. Priors for the remaining parameters are: ao ~ 
A^(0,100),/So ~ ^(0,100),/Si ~ iV(0, 100), ~ f/(0.5,10), ~ /G(2,0.1) where 

ao = E{X{s)) and /3o + /SiX(s) = iS(F(s)|X(s)). Summaries of the posterior parameter 
samples are provided in Table [TJ 


Parameter 

0.025 

Mean 

0.975 

Truth 

ao 

-0.7126 

-0.0614 

0.6214 

0 


0.2726 

0.8179 

1.3867 

0 

f3i 

0.4685 

0.5943 

0.7202 

0.5 


0.6725 

1.0692 

1.7390 

1 

4^x 

0.8242 

1.0230 

1.2124 

1.05 


0.4911 

0.8057 

1.3341 

1 

4>y 

0.8572 

1.0718 

1.3116 

1.05 


Table 1: Parameter estimates for the 2t(s) and F(s)|X(s) model. 

We consider a region centered at the location s* = (7.5, 6.5). The X(s) and R(s) values 
at this location are provided in Figure O (The full processes were realized on [0,10] x [0,10], 
but we only show a subregion here.) The interpolated surface and contour lines are produced 
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using the full 2000 locations, while the circles indicate the subset of 200 locations used to 
predict the gradient. 




i 


3 

2 

1 

0 

-1 

-2 

-3 


Figure 2: X(s) (left) and y(s) (right) subregions around s*, where we estimate the gradient. 


4 . 2 . Local Directional Sensitivity Process 

We are interested in the behavior of Zi )„y(s)/iJ„X(s) . We c onsider u = (1,0) and 
u = (0,1). Since D_^Y{s) = —D.^Y{s) flBaneriee et al.l. l2003bl) . any discussion of the 
behavior in u direction implies the opposite behavior is occurring in the opposite direction. 
When applied to the ratios, this means that the local directional sensitivity process will be 
equal for u and — u. 

Returning to the region in Figure [21 to visualize the local directional sensitivity process 
we draw samples at a grid of 125 locations denoted as {s);,..., s*}. We draw 2000 samples 
of (Vy(s^), . . . , Vy(s*), Vx(sl), . . . , Vx(s* ))' from the joint predictive distribution, again, 
given the htting data at 200 observed locations. For each of these samples we calculate 
DuY{s*)/DuX{s*). We summarize the central behavior of these Cauchy quantities at a 
given location using the median value of the ratios. 

For the two directions being considered, we plot the median predictive surface in Figure 
[3l Interpretation of these surfaces requires examination of the sign of the ratio as well as the 
magnitude. Magnitudes less than 1 suggest that the X(s) surface is changing more rapidly 
than the Y (s) surface; magnitudes greater than 1 suggest that the Y (s) surface is changing 
more rapidly than the W(s) surface; negative values suggest that one surface is increasing 
while the other decreases; positive values suggest that both surfaces are either increasing or 
decreasing. 

The direction u = (1, 0) points towards the east. The corresponding ratio surface is 
provided in the left hand plot of Figure [31 There is a peak in the ratio surface around 
(6.75,6.25), suggesting that both surfaces are either decreasing or increasing and that the 
R(s) surface is doing so more rapidly. Referring back to Figure [21 the contour lines indicate 
that both surfaces are increasing at that location looking east, and that the R(s) surface is 
doing so more rapidly. The direction u = (0,1) points towards the north. The corresponding 
ratio surface is provided in the right hand plot of Figure [3l There is a peak in the ratio 
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surface around (8.75, 6), suggesting that both surfaces are either decreasing or increasing 
and that the Y (s) surface is doing so more rapidly. Referring back to Figure [H the contour 
lines indicate that both surfaces are increasing at that location looking east, and that the 
Y (s) surface is doing so more rapidly. 
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-1 

-2 


Figure 3: Posterior median of DuY{s)/DuX{s) in the directions u = (1,0) (left) and u = (0,1) (right). 


4-3. Spatial Angular Discrepancy Process 

Consider again the region in Figure |2]and the grid of 125 locations denoted as {s^,..., s*}. 
Sampling gradients from the joint predictive distribution for (Vx(s*), Vy(s*))', we calculate 
the direction of maximum gradient as Vx (s:)/iivx(s:)ii and Vy(sn/||Vy(s:)|| for each 
sample gradient at each location in the hgure. Denote these angles (in radians) as 6'x(s*) 
and 9y{s*) respectively. 

We compute the discrepancy between these angles at each location as in Section 13.31 
and provide the posterior median values in Figure 01 Most of the region has an associated 
distance of 0, suggesting that both X(s) and R(s) are typically increasing most rapidly in 
the same direction. However, there are a few small regions where the distance peaks towards 
a value of 2, locations where the X(s) and F'(s) surfaces are increasing in nearly opposite 
directions. 

5. Duke Forest Point Pattern, Elevation Example 
5.1. Data 

Our illustrative data set is a collection of point patterns of tree species present at the 
Blackwood site in the Duke Forest in Durham, NC. The site is 5 hectares in area and exhibits 
a range of elevation. A road and powerline separate the site into three subregions, and we 
focus on the northwestern of these subregions. We consider two tree species. Flowering 
Dogwood (Cornus florida) and Sweetgum [Liquidambar styraciflua), with regard to their 
respective point patterns of locations within the site in the year 2000. Some trees have 
multiple stems observed at a single location; however we treat these as a single observed tree 
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Figure 4: Posterior median disc(s). 


at the given location. The point pattern for Sweetgum consists of 531 trees, and the point 
pattern for Flowering Dogwood consists of 570 trees. Elevation is recorded at each location 
where a tree of any species was observed, resulting in 5654 elevation observations. 

Figure |5] provides a heatmap of the elevation data. There is a clear increase in elevation 
across the region in a roughly southeastern direction. Figure [5] also provides the observed 
point patterns for each of the species. Flowering Dogwood is well dispersed across the entire 
region while Sweetgum is more abundant in the northwestern half of the region. 
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Figure 5: From left to right, observed point patterns for Flowering Dogwood and Sweetgum and observed 
elevation. 


5.2. Point Pattern Model given Elevation 

Given the observed elevation, X(s), we model the intensity for each species as a log- 
Gaussian Gox process: A(s) = exp (/So + /^i-^(s))Ao(s), Ao(s) = exp(w;j(s)), and tc^(s) ~ 
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G'F(0, p(-|02, cr^)). Placing priors on the parameters, we have the following model: 


[Y\l3o,Pi,{Xo{s),s e D}][Xo{s)\a^,(j)][Po,/3i][(j)][a^] 
= [^1/^0, A,{w^^(s),s G D}][w^{s)\a^,(j)][^o,^i][(j)][a^] 


where Y is the observed point pattern and Ao(s) is eqnivalently considered using ic^(s). 

We approximate the likelihood by dividing the region into a hne grid with cells {Ap I = 
1,..., L}. This gives us the likelihood 

L(A(s), se D] si, S2, ..., s„) Ri niA(si) exp(-A(F>)) 

XiD)^j:iexp{X'{Ai)f3 + wMi)) 


where Wz{Ai) corresponds to a realization from a Gaussian process evaluated at a representa¬ 
tive point in ea ch grid cell Ai . This likeli hood can be sampled using elliptical slice sampling, 
as described in iMurrav et ahl (120091) and iMurrav and AdamsI fl20101 1. 


5.3. Conditional Bivariate Model 

In terms of Z{s) = log(A(s)), we immediately have a conditional bivariate GP model 
with elevation: Z(s)|X(s) = /So + /SiX(s) -1- tc^(s) and X(s) = + tCx(s) with tc^(s) ~ 

GF(0, (f)z)) and Wx{s) ~ GF(0, p(-|cr^, 0a;)). Hence, we are in the framework developed 

above and can apply the proposed gradient analyses. The difference in this case will be that 
the Gaussian response Z{s) is latent and thus unobserved. The uncertainty about Z{s) is 
propagated through the model by drawing a posterior sample of the Z{s) surface for each 
posterior sample of the parameters. 
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Figure 6: Posterior median of the intensity surface for Flowering Dogwood (left) and Sweetgum (right). 


Table [2] provides the htted parameter values for elevation. Tables [3] and 0] provide the 
htted parameter values for each of the speci es. Note tha. t the (p z parameter is hxed at the 
minimum contrast estimate, as suggested in iMoller et ahl (1199811 . to facilitate identihability 


17 






















in the fitting of the log-Gaussian Cox process. The hxed values are 02 = 0.1063 and 0.0434 
for Flowering Dogwood and Sweetgum respectively. The 95% credible interval for /3i contains 
zero for Flowering Dogwood, but not for Sweetgum. This negative coefficient suggests that 
the intensity of Sweetgum decreases as elevation increases, while the intensity of Flowering 
Dogwood is not responsive to elevation changes at this scale. Figure |6] provides the posterior 
median intensity for each of the species. Both intensities have fairly low values across most 
of the domain, with a few regions of higher intensity. 


Parameter 

0.025 

Mean 

0.975 

cto 

166.1043 

167.8577 

169.4328 


7.1681 

9.4117 

12.9494 

4^x 

0.0791 

0.0892 

0.0979 


Table 2: Parameter estimates for elevation. 


5.4- Local Directional Sensitivity Process 


Parameter 

0.025 

Mean 

0.975 

Po 

-3.6416 

-3.5478 

-3.4590 

/3i 

-0.0866 

-0.0284 

0.0348 


0.3367 

0.5237 

0.7648 


Table 3: Parameter estimates for Z(s)|jy(s) model for Flowering Dogwood. 


Parameter 

0.025 

Mean 

0.975 

Po 

-4.2161 

-4.0297 

-3.8728 

Pi 

-0.4094 

-0.2619 

-0.1024 


0.7942 

1.2652 

1.9160 


Table 4: Parameter estimates for Z(s)|X(s) model for Sweetgum. 

The methods developed in previous sections allow for straight forward examination of 
DulogA(s)/DuX(s) = DuZ(s)/DuX(s), although interest is more likely in the behavior of 
the intensity surface itself. Applying the spatial gradient chain rule discussed in Section 13.41 
gives us DuA(s) = exp(Z(s))DuZ(s), with the directional derivative ratio iAuA(s)/iAuX(s) = 
exp(Z(s))DuZ(s)/iAuX(s). As before, we can simplify this in terms of the independent 
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Gaussian processes: 


DuXjs) 

DuX{s) 


,ry, ^^^U^(S) 

exp{Z{s))[PiD^Wa;{s) + D^w^{s)] 


exp(Z(s))[/? 


DuWx (s) 

DuWzjs) 

DuWxis) 


As before, (3i will center the Cauchy random variable, but now there will also be scaling 
according to the value of exp(Z(s)). 

In Figure |5] we saw a clear increase in elevation in a roughly southeastern direction. We 
approximate this direction by the unit vector u = (0.8508, —0.5255). Again, via the chain 
rule, we can consider the behavior of the directional derivative ratios in this direction for 
each of the species. Figure [7] plots the resulting posterior median Zi)uA(s)/Zi)uX(s) surfaces. 

For Flowering Dogwood the majority of the domain has a ratio close to zero. This 
suggests that the changes in the intensity are negligible compared to the changes in elevation. 
Recalling the fairly even spread of the trees in the region, as well as the non-signihcant f3i, 
this pattern makes sense. 

For Sweetgum, virtually the entire region has a negative directional derivative ratio. 
This aligns with our interpretation of the signihcantly negative (3i coefficient, namely that 
as elevation increases the intensity decreases. There are a few subregions where the change 
in intensity occurs more rapidly than elsewhere, and there is a larger subregion where the 
change in intensity is zero due to an absence of trees. The cause for subregions of rapid 
change could be further illuminated through examination of other factors in those regions. 
Similarly, the region of zero change in intensity could be roughly interpreted as Sweetgum 
having an aversion to elevations beyond a certain value. 


5.5. Spatial Angular Discrepancy Process 

Finally, we can compare the intensity and elevation surfaces by computing the discrepancy 
between their directions of maximum gradient at each location, i.e., the posterior median of 
disc{s) = 1 — cos(6'x(s) — 6*y(s)) across the region. Values close to 2 suggest the surfaces 
are most rapidly increasing in opposite directions; values close to 0 suggest the surfaces are 
most rapidly increasing in the same direction. The posterior median discrepancy surfaces 
are provided in Figure [H] for Flowering Dogwood and Sweetgum. 

The discrepancies for Flowering Dogwood roughly range between 0.8 and 1.5. There 
is no clear pattern, which supports there being no strong relationship between Flowering 
Dogwood intensity and elevation. 

The pattern for Sweetgum is quite different. All of the discrepancies appear to be between 
1.5 and 2, with most around 1.9. This suggests that the Sweetgum intensity and elevation 
are increasing in nearly opposite directions virtually everywhere in the domain. This again 
confirms the negative relationship, and additionally highlights this pattern as being slightly 
weaker in the northern part of the region. 
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Figure 7: Posterior median of DuX{s)/D^X(s) for Flowering Dogwood (left) and Sweetgum (right); u = 
(0.8508,-0.5255). 

6. Summary and Conclusions 

We have developed methodology for performing spatial sensitivity analysis for a bivariate 
process where one variable is treated as a response to the other variable. Consideration of the 
associated directional derivatives can be done jointly and results in a multivariate Gaussian 
process directly derivable from the model for the parent process £t in a Bayesian framework. 
Utilizing the posterior draws of the process parameters, all gradient analysis occurs post 
model htting. 

Using the directional derivatives, we proposed two derived processes in order to learn 
about the relationship between the response and covariate processes. The hrst is the local 
directional sensitivity process, inspired by quantities explored in standard sensitivity analy¬ 
sis. This process captures local variation in the relationship between the two variables and 
provides deeper insight into their relative behavior across the region. The second is the spa¬ 
tial angular discrepancy process, capturing the discrepancy between the directions in which 
the process surfaces are most rapidly increasing. Spatial plots of this discrepancy surface 
highlight regions of the domain where the two processes behave most similarly and most 
differently, again informing on the variables’ spatial relationship. 

Our application involved elevation data and point patterns of trees collected from the 
Duke Forest. Using a log-Gaussian Cox process, we studied local directional sensitivity of 
the intensity to elevation through the use of a spatial gradient chain rule. Through the two 
tree species, this example illustrates the different results one would expect from a spatial 
sensitivity analysis when the variables are signihcantly related versus when they are not. 

The current theory provides opportunity for several extensions and applications. Many 
ecological data sets are observed at multiple time points, in part to see if the relationships 
between the variables of interest are changing over time. With this in mind, future work 
on gradient analyses may involve the incorporation of temporal effects. Ecological data sets 
can also have multiple responses such as leaf traits that are being related to multiple climate 
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Figure 8: Posterior median disc{s) for Flowering Dogwood (left) and Sweetgum (right). 


covariates (e.g. iThuiller et a1 


2nn4 . any or all of which may have relationships that could 
be better highlighted through a spatial gradient analysis under a joint model. Finally, the 
spatial gradient chain rule will allow us to consider novel non-Gaussian responses modeled 
through latent Gaussian process models as described in Section [331 


Appendix A. 

We show that D^Wy^s)/D^w^i^s) is a well defined spatial stochastic process. 

First consider two locations s and s'. For notational convenience, let DuWy{s) = rii, 
DuWy{s') = 77-2, DuWx{s) = mi, and DuWx{s') = m 2 . We can write the joint distribution of 
the ratios in terms of the joint distributions for the two Gaussian processes: 


P{— < Ti, —^ < r2) = P{ni < rimi, n2 < r2m2, mi > 0, m2 > 0) 

mi m 2 

+ P{ni < Timi, 712 > ’" 2 ^ 2 , 'mi > 0, m 2 < 0) 

+ P{ni > riTTii ,772 < 'r2'm2, mi < 0, m 2 > 0) 

+ P(77i > riTTli, 772 > ^"2^25 < 0, 7772 < 0) 


f*oo poo primi pr2m2 


fO Jo 

poo pO 


r*rimi poo 


+ 


/i^(7ii, 712 )/g(?^i, m2)dn2dnidm2dmi 


/ii-(77i, n2) fci^i, m2)dn2dnidm2dmi 


' —00 J —00 
noo /‘CO 


' r2m2 
(*r2m2 


‘ —CO Jo J rimi J —CO 


fO 


^•0 


/ii-(77i, n2) fci^i, m2)dn2dnidm2dmi 


/i^(77i, n2)fG{mi,m2)dn2dnidm2dmi 


' —CO J —CO J rimi J r 2 m 2 
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The above joint distribution depends only on integrals of multivariate normal densities, so 
consistency across permutations of the labels is clearly satisfied. 

At a single location the directional derivative ratio is a ratio of two independent normal 
random variables. As such, it will have a Cauchy distribution with a scale parameter depen¬ 
dent on the parameters in the covariance functions K{-) and G{-), namely SD{ni)/SD{mi). 
Next we show that marginalizing the bivariate distribution over r 2 will reduce to the known 
distribution for the univariate case: 


P(—<ri) = 

mi 


+ 


F(— < ri, — < r 2 )dr 2 

mi m 2 


' —00 

^00 roo nOG nrimi /‘T'2m2 

00 ./O Jo J —00 J —00 

POO pO poo primi poo 


/ii-(ni, n2)fG{'mi, m2)dn2dnidm2dmidr2 


fxini^ n2) TT^2)dn2dnidm2dmidr2 


'—<X) J —00 J ^ J —00 Jr 2 m 2 

poo pO poo poo pr 2 m 2 


+ 


+ 


' —00 J —CO ^0 J rimi J —CO 
poo /*0 /•O poo poo 


' —CO J —CO J —CO J rimi J r 2 m 2 
poo poo primi 


+ 


' 0 Jo J —CO 
poo pO primi 

0 J —CO J —CO 
/•O poo poo 


fK{ni,n2)fG{'mi,m2)dn2dnidm2dmidr2 
/x(n.i, n2)/G(’^i5 m2)dn2dnidm2dmidr2 
fK{ni)fG{mi, m2)dnidm2dmi 
/x(’^i)/g(^i, m2)dnidm2dmi 
fK{ni)fG{mi,m2)dnidm2dmi 


+ 


' —CO Jo J rimi 

pO rO poo 

' —CO J —CO Jr\m\ 

poo primi 


fK{ni)fG{mi, m2)dnidm2dmi 

pO poo 


fkini)fGimi)dnidmi + 


fk(,ni)fG{mi)dnidmi 


' 0 J —CO 


' —CO J rimi 


We can then rewrite this as: 

Hi 

P {— < ri) = / F^(rimi)/G(mi)dmi 

mi Jo 

poo 

FKinmi) fG{mi)dmi 


FK{-rimi)fG{mi)dmi 


= 2 


The a ssociated density will then be 2 mifK{rimi)fG{mi)dmi, which was shown in [A ndrews and Mallow: 
(119741 ) to be a Cauchy distribution. In this case we will have a Cauchy with scale parameter 
SD{ni)/SD{mi), as desired. 

The marginalization is straight forward for larger dimensions with marginalization over 
the two normal densities occurring in a similar way. 
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Appendix B. 

We are interested in the marginal behavior of f{0x (s), 6*y( s)) at a location s where 
tan(6'x(s)) = L>(o,i)A(s)/L)(i,o)^(s), tan(6'y(s)) = L)(o,i)F(s)/L)(i,o)>^(s), and 
(Vx(s), Vy(s))' = (L)(i,o)y(s),D(o,i)y(s),D(i,o)A(s),L)(o,i)X(s))' are normally distributed 
with the covariance structure provided in Equation 1. 

After converting to polar coordinates, the change of variables formula gives the following 
integral, where g{-) is the multivariate normal density associated with (Vy (s), Vx(s))^ For 
clarity we suppress the index (s): 


POO poo 


f{Ox,0Y) = 


'0 ^0 

roo 


g{ry cos 9y, Vy sin ^y, cos 6'x, sin 6x)ryrxdrydr\ 


exp(--(ac- A )r^) 


lo "^(2^ 


X 


ry-j=exp{--{ry - (dr^ cos Oy)dry ) dr^ 


1 

roo ^ 1 

^/WF\ 

|S| J 

0 \/^ \ 


+ 


Tx exp(--(ac - Ayrydr, 


1 A 


\/ (27r)2|S| ^ \/2'Ka \fd Jq 


r^^iAryV:^ exp{--{ac - AAr^)dr, 


+ 


1 ~ 1 ~ 

/ — ^(j){Aryexp{--{ac-A‘^)ri)dr, 

lo ayzTi ^ 


1 


-y/(27r)2|S| y av^ \ ac(ac-d2) {ac-A^)^/'^ 




+ 


m 


1 ( A ( </>(0)A , V27r K _ T (0 f) a / \ 

Y^( 27 r )2 |S| I \.ac{ac—A‘^) {ac—A‘^)^/‘^ ^ ^ ’ y ac ^ ' J a^c\/^ 


yJ‘l'K 


’ '/ ac ' 




A > 0 
, A < 0 


where |S| = {alcplY{a^cplf, a = c = {a^c^l + f3'^c^lal)/{alc^l), A = ^/d cos(9x - 

dy), and L(0,0,p) is the zero mean bivariate normal cdf with correlation p and standard 
deviations equal to 1 evaluated at (0, 0)'. 

Several integration steps were treated as the expected value of a truncated normal distri¬ 
bution. An additional integration step required integration by three parts, then computation 
of the probability that the sum of a normal and a truncated normal are less than some con¬ 
stant. This probability was available in closed form in the query by Lipow et ah f 19641) and 
is written in terms of the bivariate normal cdf. Somewhat simplihed, this gives the density 
provided in the text: 


f{ex,ey) 




y/27vA 


a(27r)3/2y^ \^ac(ac-A2) (ac-A2)3/2 

_ A^(/)(0) I \/2^A 

a(2'K 




Ai'i I Ao) '\ 


t)3/2W[s| \ ac{ac—A^) [ac—Ap!^ \ ' ’ 


I ) 


'hi') ] + ^ 

ac ' j ac 


A > 0 
A < 0 
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